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Abstract 

We discuss a recently proposed analytic solution to the Thomas- 
Fermi (TF) equation and show that earlier approaches provide more 
accurate results. In particular, we show that a simple and straightfor- 
ward rational approximation to the TF equation yields the slope at 
origin with unprecedented accuracy, as well as remarkable values of 
the TF function and its first derivative for other coordinate values. 



1 Introduction 

The Thomas-Fermi (TF) equation has proved useful for the treatment of 
many physical phenomena that include atoms [1—5], molecules [3,6], atoms in 
strong magnetic fields [1,4,5], crystals [7] and dense plasmas [8] among others. 
For that reason there has been great interest in the accurate solution of that 
equation, and, in particular, in the accurate calculation of the slope at origin 
[9-11]. Besides, the mathematical aspects of the TF equation have been 
studied in detail [14,15]. Some time ago Liao [16] proposed the application 
of a technique called homotopy analysis method (HAM) to the solution of 
the TF equation and stated that "it is the first time such an elegant and 
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explicit analytic solution of the Thomas-Fermi equation is given". This 
claim is surprising because at first sight earlier analytical approaches are 
apparently simpler and seem to have produced much more accurate results 
[10-13]. Recently, Khan and Xu [17] improved Liao's HAM by the addition 
of adjustable parameters that improve the convergence of the perturbation 
series. 

The purpose of this paper is to compare the improved HAM with a 
straightforward analytical procedure based on Pade approximants [13] sup- 
plemented with a method developed some time ago [18-26]. In Section [2] 
we outline the main ideas of the HAM, in Section [3] apply the Hankel-Pade 
method (HPM) to the TF equation, and in Section H] we compare the HAM 
with the HPM and with other approaches. 



2 The homotopy analysis method 

In order to facilitate later discussion we outline the main ideas behind the 
application of the HAM to the TF equation. The TF equation 



u "{x) = y^p> "(0) = 1, w(oo) = (1) 

is an example of two-point nonlinear boundary-value problem. When solving 
this ordinary differential equation one faces problem of the accurate calcula- 
tion of the slope at origin u'(0) that is consistent with the physical boundary 
conditions indicated in equation (flj) . 

In what follows we choose the notation of Khan and Xu [17] whose ap- 
proach is more general than the one proposed earlier by Liao [16]. They 
define the new solution g(£) = ju(x), where £ = 1 + Ax and rewrite the TF 
equation as 

(£ - 1)A V(0 3 - <?(0 3 = (2) 

where 7 is the inverse of the slope at origin (V(0) = I/7) and A is an 
adjustable parameter. Khan and Xu [17] state that the solution to Eq. §Z§ 
can be written in the form 

00 

9(0 X-V ' (3) 
that reduces to Liao's expansion [17] when A = 1. 
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In principle there is no reason to assume that the series converges and 
no proof is given in that sense [16, 17]. Besides, the partial sums of the series 
([3]) will not give the correct asymptotic behaviour at infinity [14,15,27] as 
other expansions do [10, 11]. 

Liao [16] and Kahn and Xu [17] do not use the ansatz ([3]) directly to 
solve the problem but resort to perturbation theory. For example, Kahn and 
Xu [17] base their approach on the modified equation 

(i - q )c im q) - 9o(0) = qtM 9), r( g )] (4) 

where £ and M are linear and nonlinear operators, respectively, < q < 1 
is a perturbation parameter and h is another adjustable parameter. Besides, 
<7o(0 is a conveniently chosen initial function and $(£; q) becomes the solu- 
tion to equation when q — 1 [17]. Both $(£; q) and T(q) are expanded in a 
Taylor series about q = as in standard perturbation theory, and T(0) = 70 
is another adjustable parameter [17]. 

The authors state that HAM is a very flexible approach that enables one 
to choose the linear operator and the initial solution freely [16, 17] and also 
to introduce several adjustable parameters [17]. However, one is surprised 
that with so many adjustable parameters the results are far from impressive, 
even at remarkable great perturbation orders [16,17]. For example the [30/30] 
Pade approximant of the HAM series yields u'(0) with three exact digits [17], 
while the [1/1] Pade approximant of the 5 expansion [28] provides slightly 
better results [29,30]. A more convenient expansion of the solution of the 
TF equation leads to many more accurate digits [10,11] with less terms. 



3 The Hankel— Pade method 



In what follows we outline a simple, straightforward analytical method for 
the accurate calculation of u'(0). In order to facilitate the application of the 
HPM we define the variables t = x 1/2 and f(t) = u{t 2 f/ 2 , so that the TF 
equation becomes 



T(f,t) = t f(t)f'(t) + f(ty -f( t )f(t)-2t 2 f(ty 







(5) 



We expand the solution f(t) to this differential equation in a Taylor series 
about t = 0: 

/(*) = E // (6) 
3=0 
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where the coefficients fj depend on / 2 = /"(0)/2 = u'(0)/2. On substitution 
of the series (jSD into equation <^ we easily calculate as many coefficients fj 
as desired; for example, the first of them are 

fo — 1) /i = 0, fz = -, A = — /s = — • • ■ (7) 

The HPM is based on the transformation of the power series into a 
rational function or Pade approximant 

WW) = (8) 

One would expect that M < N in order to have the correct limit at infinity; 
however, in order to obtain an accurate value of / 2 it is more convenient to 
choose M = N + d, d = 0, 1, . . . as in previous applications of the approach 
to the Schrodinger equation (in this case it was called Riccati-Pade method 
(RPM)) [18-26]. 

The rational function (jSJ) has 2N+d+l coefficients that we may choose so 
that T([M/N], t) = 0{t 2N+d+1 ) and the coefficient / 2 remains undetermined. 
If we require that T([M/N],t) = 0(t 2N+d+2 ) we have another equation from 
which we obtain / 2 . However, it is convenient to proceed in a different (and 
entirely equivalent) way and require that 

2N+d+l 

[M/N](t)~ fjt j = 0(t 2N+d+2 ) (9) 

j=0 

In order to satisfy this condition it is necessary that the Hankel determinant 
vanishes 

H D = l/i+J+d+llij=0,l,..JV = °> ( 10 ) 

where D = N + 1 is the dimension of the Hankel matrix. Each Hankel 
determinant is a polynomial function of / 2 and we expect that there is a 
sequence of roots f^'^, D = 2, 3, . . . that converges towards the actual value 
of u'(0)/2 for a given value of d. We compare sequences with different values 
of d for inner consistency (all of them should give the same limit). Notice 
that a somewhat similar idea was also proposed by Tu [12], although he did 
not develop it consistently. 

Present approach is simple and straightforward: we just obtain the Taylor 
coefficients fj from the differential equation (jSj) in terms of / 2 , derive the 
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Hankel determinant, and calculate its roots. Since is the first nonzero 
coefficient that depends on f 2 we choose Hankel sequences with d > 3. 

The Hankel determinant Hf) exhibits many roots and their number in- 
creases with D. If we compare the roots of Hf) with those of Hf ) _ l we 
easily identify the sequence f? A that converges towards the actual value 
of f 2 . Fig. [T] shows log 2/2°'^ — 2f2 D ~ 1,d ^ for D = 3, 4, . . . that provides 
a reasonable indication of the convergence of the sequence of roots. We 
clearly appreciate the great convergence rate of the sequences with d = 3 
and d — 4. For example, for d = 3 and D < 30 it is approximately given 
by 2f 2 [D ' 3] - 2/j D_1 ' 3] | = 14.2 x i(r - 705D From the sequences for D < 30 
we estimate u'(0) = —1.58807102261137531 which we believe is accurate to 
the last digit. We are not aware of a result of such accuracy in the literature 
with which we can compare our estimate. It is certainly far more accurate 
than the result obtained by Kobayashi et al [9] by numerical integration that 
is commonly chosen as a benchmark [16, 17]. 

Present rational approximation to the TF function is completely different 
from previous application of the Pade approximants, where the slope at origin 
was determined by the asymptotic behaviour of at infinity [13] . Our approach 
applies to u(x) 1 ^ 2 and the slope at origin is determined by a local condition 
at that point which results in the Hankel determinant ffTTJ|) . In this sense 
our approach is similar to (although more systematic and consistent than) 
Tu's one [12] as mentioned above. 

Once we have the slope at origin we easily obtain an analytical expression 
for u(x) in terms of the rational approximation (jSJ) to /(£). In order to have 
the correct behaviour at infinity we choose N = M + 3 [13]. Table [T] shows 
values of u(x) and its first derivative for 1 < x < 1000 (the approximation 
is obviously much better for < x < 1) given by the approximant [5/8]. 
Our results are in remarkably agreement with the numerical calculation of 
Kobayashi et al [9] and are by far much more accurate than those provided 
by the HAM [16,17]. Notice that we are comparing a [5/8] Pade approx- 
imant on the straightforward series expansion ([6]) with [50/50] and [30/30] 
approximants on an elaborated perturbation series [16, 17]. 



4 Conclusions 

Any accurate analytical expression of the solution u(x) to the TF equation 
requires an accurate value of the unknown slope at origin u'(0), and the HPM 
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provides it in a simple and straightforward way. In this sense the HPM ap- 
pears to be preferable to other accurate approaches [9-11] and is far superior 
to the HAM [16,17]. Notice for example that our estimate 2/j 5,3] = -1.588, 
based on a rational approximation [7/4], is better than the result provided by 
a [30/30] Pade approximant on the improved HAM perturbation series [17]. 
Besides, by comparing Table 2 of Khan and Xu [17] with our Fig. [1] one 
realizes the different convergence rate of both approaches. One should also 
take into account that the HPM does not have any adjustable parameter for 
tuning up its convergence properties, while, on the other hand, the "flexible" 
HAM with some such parameters plus a Pade summation results in a much 
smaller convergence rate [16,17]. 

We also constructed a Pade approximant [5/8] from the series §6§ and ob- 
tained the TF function and its derivative with an accuracy that outperforms 
the [50/50] and [30/30] Pade approximants on the HAM perturbation se- 
ries [16,17]. It is clear that the HPM is by far simpler, more straightforward, 
and much more accurate than the HAM. 

In addition to the physical utility of the HPM we think that its mathe- 
matical features are most interesting. Although we cannot provide a rigorous 
proof of the existence of a convergent sequence of roots for each nonlin- 
ear problem, or that the sequences will converge towards the correct phys- 
ical value of the unknown, a great number of successful applications to the 
Schrodinger equation [18-26] suggest that the HPM is worth further inves- 
tigation. Notice that we obtain a global property of the TF equation w'(0) 
from a local approach: the series expansion about the origin §6§. The fact 
that our original rational approximation ([8]) does not have the correct be- 
haviour at infinity is not at all a problem because we may resort to a more 
conventient expansion [13] once we have an accurate value of the unknown 
slope at origin. 

Finally, we mention that the HPM has recently proved successful for the 
treatment of other two-point nonlinear equations [31] of interest in some 
fields of physics [32-34]. 
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Table 1: Values of the Thomas- Fermi function and its derivative 



X 


u(x) 


—u'{x) 


1 


0.424008 


0.273989 


5 


0.078808 


0.023560 


10 


0.024315 


0.0046028 


20 


0.005786 


0.00064727 


30 


0.002257 


0.00018069 


40 


0.001114 


0.00006969 


50 


0.000633 


0.00003251 


60 


0.000394 


0.0000172 


70 


0.0002626 


0.000009964 


80 


0.0001838 


0.000006172 


90 


0.0001338 


0.000004029 


100 


0.0001005 


0.000002743 


1000 


0.000000137 


0.00000000040 




Figure 1: log 2ff A - 2}f~ lA for d = 3 (circles) and d = 4 (filled circles) 
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